Electronic Wavefunction in Empirical Tight-Binding Theory 
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A tight-binding procedure to calculate the single-electron wavefunctions of zinc-blende semicon- 
ductors is proposed. It is based on linear combinations of Slater-type orbitals whose screening coef- 
ficients are self-consistenly extracted from the optical matrix elements of the tight-binding Hamil- 
tonian. The Bloch functions obtained in the extended-basis spds* tight-binding model demonstrate 
very good agreement with first-principles wavefunctions. We apply this method to the calculation 
of the excitonic fine structure of bulk GaAs. Beyond semiconductor nanostructures, this work is 
a fundamental step toward modeling many-body effects from postprocessing wavefunctions within 
the Slater and Koster theory. 

PACS numbers: 



Electronic excitations in nanoscale structures stand at 
the heart of modern band-structure theory. In this con- 
text, several theoretical approaches have been used for 
addressing the electron-correlation problem through the 
calculation of the single-particle energies and wavefunc- 
tions of electrons and holes pp. "Screened-hybrid func- 
tionals" with GW method now allow first-principles cal- 
culations to model succesfully the bandgaps and effective 
masses of semiconductors [2 . However, they remain lim- 
ited to bulk materials or to rather small clusters. The 
scientific and technological importance of representing 
accurately the excited-state properties of semiconduc- 
tors and their nanostructures has led to the development 
of various empirical approaches, ranging from empirical 
pseudopotentials over k.p models to tight-binding (TB) 
methods. The two first ones have been extensively used 
for calculating the excitonic fine structure of nanostruc- 
tures in recent years. For example, the electron-hole ex- 
change interactions in semiconductor quantum dots has 
been thoroughly studied using atomistic pseudopoten- 
tials [3 J. This approach represents the state-of-the-art 
of nanostructure simulations but, in present implemen- 
tations, suffers from a parametric poorness that does not 
allow for an accurate description of the full Brillouin 
zone. On the other hand, tight-binding has existed for 
many years as a convenient and conceptually transpar- 
ent model for the description of electronic structure of 
molecules and solids. It provided the basis for construc- 
tion of many body theories such as the Hubbard model 
[4] and the Anderson impurity model [5 j. Slater and 
Koster (SK) called it the tight binding or Bloch method 
and their 1954 seminal paper [6] provided the system- 
atic procedure for formulating a tight-binding Hamilto- 
nian. In this method, the crystal potential is approx- 
imated as a sum of spherically symmetrical potentials 
around each atom. This allows the electronic wavefunc- 
tions to be developed on a set of atomic-like orbitals 



(called the Lowdin orbitals) having well defined angu- 
lar properties, but unknown radial dependencies. The 
Hamiltonian matrix elements between Lowdin orbitals 
are treated as "disposable constants" with which one 
can fit band structures that have been experimentally 
determined or calculated using some more accurate tech- 
niques. The fit is performed in k-space, removing any ne- 
cessity to further characterize the local wavefunctions in 
real space. Very interestingly, interaction with the elec- 
tromagnetic field is built-in without any ambiguity from 
a point of view of the gauge invariance using a derivation 
of the Hamiltonian matrix elements in momentum space: 
Pim,i'm' = iH^im I VkH \ ^i>m>) [3 E] • Optical prop- 
erties can be consequently calculated in semi-empirical 
tight-binding models without adding parameters. In this 
context, the extended-basis TB model is known to give 
a description of dielectric properties equivalent to best 
ab initio calculations within the one-electron approxima- 
tion. However, the question of on-site transition matrix 
elements, not included in this approach, is still debated 
within the TB community [7HT0]. Altogether, the model 
major qualities are the transferability of parameters from 
bulk materials to nanostructures, the unique ability to 
describe with the same accuracy electronic properties in 
any region of the Brillouin zone and the capacity to han- 
dle large supercells, up to million atoms. However, when 
it comes to interactions (in particular, short range inter- 
actions) between quasi-particles, the TB method is ham- 
pered by the lack of knowledge of the local wavefunction. 
Existing calculations of Coulomb matrix elements in crys- 
tal semiconductors [11] use atomic approximations on the 
radial dependence of the basis orbitals [12] that give a 
poor representation of dipole matrix elements. In this 
letter we solve this theoretical issue that has remained 
open ever since the seminal work of Slater and Koster, 
by defining formally the radial wavefunctions out of the 
electronic Hamiltonian. 
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We start with a set B of spds* basis functions in the 
form of normalized Slater- type orbitals (STO) & n im(r) = 
v /(2a) 2n + 1 /(2n)! Yj, m (<9, ^)r n - 1 e- QLr , where n is the first 
quantum number and a is a screening parameter [T2] . 
STOs are largely employed in quantum chemistry but do 
not fulfill the orthogonality condition, since finite overlap 
exists between two STOs localized at different sites of the 
crystal. 

Firstly, the orbital overlap matrix S is calculated 
including all orbitals up to a cut-off distance Rq. Ro is 
related to the smallest a, and is taken large enough so 
that overlap with remote atoms be negligible: This step 
is done analytically [13], and in practice, we found that 
for a = 0.5, overlap with neighbors located farther than 
3 lattice parameters (17A) can safely be neglected. Note 
that, unlike real atomic orbitals, s and s* STOs cannot 
be orthogonal on-site. The most physical solution to 
this difficulty consists in substituting s* with a Gramm- 
Schmitt combination s* = (s* - (s\s*)s)/y/l - (s^*) 2 . 
Then, an orthogonal basis B ort h can be obtained using 
the Lowdin orthogonalization procedure B ort h — 
[T4] . The orthogonalized STOs will serve as trial func- 
tions for the unknown Lowdin orbitals. The expansion of 
electronic eigenfunctions (Bloch functions) in the basis 
B is obtained by multiplying the eigenvectors of the 
Hamiltonian matrix of the sp 3 d 5 s* model by the matrix 
£2, which definitely provides their representation in 
real space. However, the 9150 x 9150 S'-matrix obtained 
this way contains considerably redundant information 
and its large size makes its inversion computation- 
ally difficult. Since in fine we are interested in the 
crystal eigenstates (Bloch functions), instead of the 
procedure sketched above, we equivalently construct, 
orthogonalize and invert a 40 x 40 S overlap matrix 
between Bloch sums of STOs, truncated to the cut-off 
distance Ro. Although this may not be obvious at first 
glance, S and S overlap matrices produce equivalent 
results. Then the momentum matrix elements are 
calculated from the Bloch functions by the relation 
Pim,i'm' = iH^im I Vr | *i'm')- Tni s derivation in real 
space involves a sum of matrix elements between two 
STOs that can be calculated analytically. We note that 
the orthogonalization procedure induces intra-atomic 
matrix elements of the momentum operator, compatible 
with the gauge invariance as discussed by Sandu [T5] . 
This, however, does not solve completely the issue of 
on-site matrix elements. Finally, the screening param- 
eters are fitted into a genetic algorithm to reproduce 
the optical matrix elements derived from the electronic 
Hamiltonian in k-space for consistency. The consider- 
ation of different bands and different high symmetry 
points in the Brillouin zone provides more than necessary 
information for the fit convergence. However, we observe 
that, despite the efficiency of the genetic algorithm, 
rather different sets of screening parameters can give 
similar "fitness" parameters. This difficulty is mostly 



TABLE I: Optimized Slater orbital screening coefficients for 
gallium (Ga) and arsenic (As), compared with Slater's atomic 
screening coefficients [12] 



Orbital 




Ga 




As 




Ref [I2] 


This work 


Ref [H] 


This work 


4s 


1.35 


1.83 


1.7 


1.94 


4p 


1.35 


1.77 


1.7 


1.79 


4d 


0.27 


0.93 


0.27 


0.96 


5s 


0.32 


1.64 


0.4 


1.74 



related to the fact that the dispersion of the upper bands 
cannot be exact, due to the non-completeness of finite 
basis. It is cured by constraining the parameter space 
in such a way that the d(T 12) orbital, which is close to 
a free-electron state, agrees with independent empirical 
pseudo-potential or ab initio calculations. 

We applied the procedure explained above to the pro- 
totype systems of Ge and GaAs. The electron configu- 
ration of Ge is [Ar]3d 10 4s 2 4p 2 . In the spds* model, the 
deep 3d states are discarded and the basis is formed by 
the orbitals 45, the three orbitals 4p, the five empty or- 
bitals 4d and the empty orbitals 5s. When building the 
STO basis S, we keep fixed the first quantum number n 
of these orbitals and introduce one adjustable screening 
parameter a for each symmetry type. Alternatively, as 
often done in quantum chemistry, we can improve para- 
metrical flexibility by considering that each element of 
the starting basis is a linear combination of q STOs in- 
stead of one. This does not change much the model, 
but increases to Aq the number of fitted parameters. For 
GaAs, since there are two different atoms in the unit cell, 
the number of parameters is twice that for Ge. The fit- 
ted screening parameter for Ga and As are given in Table 
[TJ and contrasted with the Slater atomic screening con- 
stants. Those for Ge are close to averaged values for Ga 
and As. At the end of the fitting procedure, the global 
discrepancy on the sum of all interband matrix elements, 
calculated at the T, X and L points of the Brillouin zone, 
is less than 15% with one Slater orbital per atomic state 
and less than 7% with a linear combination of two Slater 
orbitals for each atomic state. By changing the relative 
weights of different spectral or Brillouin zone regions in 
the genetic algorithm fitness function, the discrepancy 
at e.g. T can be minimized down to the percent range. 
The residual discrepancy has three distinct physical ori- 
gins: i) the difference between orthogonalized STOs and 
the actual Lowdin orbitals, which can, to some extent, 
be minimized by increasing the number of STO compo- 
nents per orbital in the basis, ii) the lack of complete- 
ness of the spds* basis, and hi) the missing intra-atomic 
contribution in the k-space derivation method. Table [TT] 
shows the main momentum elements Po = — i(s c \p x \x v ) , 
Pi = — i(s c \p x \x c ) and Qo = —i{x c \p y \z v ) obtained for 



TABLE II: Main interband matrix elements (in eV A) for Ge 
and GaAs at the T point, calculated from differents models 







WFl a 


WF2 6 


LDA+GW C 


Hamiltonian^ 


Ge 


Po 
Qo 


7.69 
8.29 


10.18 
8.42 


8.49 
7.32 


10.14 
8.70 


GaAs 


Po 
Pi 
Qo 


7.38 
0.80 
8.16 


9.88 
0.93 
8.31 


8.35 
1.38 
7.37 


9.82 
0.11 
8.72 



ab real space calculation from TB wavefunctions with respectively 
one and two Slater orbital for each basis element 

c real space calculation from LDA+GW wavefunctions, ABINIT 
code 

^calculated from Hamiltonian derivation 



Ge and GaAs at the T point. For P and Q, all results 
compare well with those of the 14-band k.p model [TB]. 
Note that for GaAs the Slater's atomic screening coef- 
ficients give: Po = 1-5 eVA and Qo = 5.3 eVA, which 
underlines a poor description of the local wavefunctions 
in real space. A discrepancy is observed for Pi between 
the real space calculations and Hamiltonian derivation 
[T7] . This might be a trace of the methodological limita- 
tion relative to intra-atomic contributions. While further 
work is required to explore the method limitations and 
improve the results, the present achievement is already 
sufficient for most practical purposes. 

Once screening parameters best reproducing the inter- 
band matrix elements are obtained, the different Bloch 
functions can be plotted and compared with ab initio cal- 
culations. First-principle calculations of single-electron 
wavefunctions are performed with ABINIT code [TSJ HH] 
in the local density approximation (LDA) for exchange 
and correlation, completed by a self-consistent GW cor- 
rection. A plane-wave set with a wave vector cutoff of 
9 a.u. is used to expand the wavefunctions and the recip- 
rocal space integration is done over 8x8x8 Monkhorst- 
Pack grid. Figure [T] shows valence band states s v and 
y' v = y v — x v , in both TB and ab initio calculations. 
The overall quantitative agreement is very good, since the 
overlap between TB and ab initio densities is always bet- 
ter than 95%. Yet, TB wavefunctions appear somewhat 
less localized in the sense that they have larger density in 
regions where the ab initio density is almost zero, and the 
Ga / As asymmetry is more pronounced in the ABINIT 
result. The most significant difference is for the deep s v 
state near the atomic sites, for which TB density is sig- 
nificantly smaller. This probably reflects the difference of 
projection basis between the two models: TB wavefunc- 
tions are expanded in the basis of Slater orbitals which 
have a node at the atomic sites while the ABINIT wave- 
functions are expanded in a basis of plane- wave functions 
that may be maximum on the atomic sites. The wavevec- 
tor cut-off used in the ABINIT calculations is important, 
because this approach cannot describe the region located 
less than l/fc C ut-off from atomic sites. Yet, cut-off does 
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FIG. 1: (Color online) Isodensity contours of the s v and 
y' v = x v — y v valence Bloch function at the zone center in bulk 
GaAs in the plane (110). TB calculation (left) is compared 
with ABINIT calculations (right). 



not suffice to explain the observed difference. In our TB 
approach, finite on-site value for the s v state results from 
the contribution of neighboring atoms and it is quite sen- 
sitive to the STO screening parameters. 

In Figure [2] we show the conduction Bloch functions 
also calculated with the same two models. Again TB 
wavefunctions are very similar to those calculated in the 
LDA+GW approximation, except for a significant differ- 
ence for the s c state density in the vicinity of atomic sites. 
In order to clarify this issue, we used the SIESTA code, 
which is based on DFT expanded in a strictly localized 
orbitals set. SIESTA results for s v and s c actually agree 
very well with our TB results. We note that electron hy- 
perfine interaction constants, that are well documented, 
scale as s c 2 (r = 0) and could serve as a quantitative test. 
A most striking result is the TB ability to reproduce the 
wavefunctions of the nearly free electron states s* and d. 

The major interest of having a real space represen- 
tation of wavefunctions is the ability to study many- 
body problems. In the following, we illustrate the po- 
tential of our wavefunction derivation method by cal- 
culating the exciton binding energy and its fine struc- 
ture due to electron-hole exchange interaction. We first 
reintroduce spin-orbit interaction and make the classical 
Lowdin partition of valence band into T 8 and T 7 states. 
Coulomb interaction can be specified by the matrix el- 
ements (ra', k^; n', kjjf/^lra, k e ; n, k^), where U eh = 
e 2 /ft|r e — Th\. This Coulomb matrix element involves 
both direct and exchange terms. Here, |ra, k e ; n, k^) 
is the two-particule excited state, k is the permittiv- 
ity, and ^ m ,k e = (r|ra,k e ) and ^ n>k/i = (r\n,k h ) 
are the Bloch wave functions in electron and hole rep- 
resentations, respectively [20j [21]. The present ex- 
pansion of Bloch functions as linear combinations of 
Slater orbitals allows to expand the electron-hole inter- 
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FIG. 2: (Color online) Isodensity contours in the (110) plane 
for the s c , y' c and d(T\2) conduction Bloch functions at the 
zone center in bulk GaAs TB calculation (left) and ABINIT 
calculations (right) 

action in terms of Coulomb matrix elements between 
STOs V/ imi> / 2m2j /3 m3j i 4m4 = (3>/ imi (r — Ri), 3>/ 2m2 (r — 
R2)|^ e/l |^ 3 m 3 (r - R 3 ), $i 4 m 4 (r - R 4 ))- Restricting the 
expansion to two-center contributions (Ri = R3 and 
R2 = R4), the evaluation of the integrals can be done 
quasi analytically using the expansion of the Coulomb 
potential in terms of spherical harmonics centered on the 
same site when Ri = R 2 [22 J, and a bipolar expansion 
when Ri ^ R 2 [23 . Following [H [25], we introduce 
a r-dependence of the permittivity keeping unscreened 
the on-site and first neighbors interactions while screen- 
ing interactions between distant atoms with the mate- 
rial static dielectric constant. Then, we solve the Bethe- 
Salpeter equation (BSE) [26-32 , expressed in terms of 
calculated electron-hole interactions and single particle 
energies. The BSE is an eigenvalue problem of infinite 
dimensionality: 

(^c,k+Q/2 - ^,k-Q/2)^ck 

+ Sv BZ d * k, Y, v >,Av^u* h \v f c'U)A v , c/ u = n s A vcli 

Where ^ c ,k+Q/2 an d ^,k-Q/2 are the electron and 
hole energies respectively. Resolution of BSE gives 
the exciton wavefunction components A vc \^ and the 
excitation energies Qs- To make the problem tractable 
continuous integration with respect to k' was replaced 
by a discrete scheme. Following [33j [34], to calculate 
the exciton spectra and binding energie at the T point, 
the integration was performed over a small region near 




Q// [111] (27t/a) 

FIG. 3: (Color online) Dispersion of the exciton states of 
bulk GaAs for Q along the [111] direction. The doted line 
at small Q is an extrapolation of calculated points: the long 
range exchange is exactly zero at Q = 0, and builds up for 
very small values of Q , for which convergency is more difficult 
to ensure. The wavevector of light, for which strong polariton 
features would add to the present picture, is indicated with 
a vertical line. The inset shows schematically the different 
contributions to the fine structure of Tgv x T6 C fundamental 
exciton. 

the position of the band extrema (|k| < 0.015 2tt/o). 
This region was divided into a 11 x 11 x 11 uniform 
grid. For an exciton wave vector Q = 0, we find an 
excitonic binding energy — 4.75 meV. In addition, 
the eightfold degenerate T 8v x r 6c fundamental excitonic 
transition is split by short range exchange interaction 
into one twofold, and two threefold degenerate excitons. 
The twofold and threefold J = 2 "dark excitons" are 
split by ^ an i s = 0.02 \ieV . This anisotropy splitting 
is due to the zinc-blend structure which does not 
allow more than threefold degeneracy. We expectedly 
find a very small value for £ a nis- The J = 1 "bright 
exciton" threefold state is separated from the J = 2 
states by the short range exchange splitting A exc . We 
get Aexc — 20.6 /ieV, in agreement with experimental 
determination [35]. When one moves away slightly 
from Q = 0, the J = 1 excitons are further split 
by the long range exchange interaction into twofold 
degenerate, optically active transverse excitons, and 
a longitudinal exciton. The energy difference corre- 
sponds to the longitudinal-transverse splitting Alt, for 
which we find a value Alt = 105.3 fieV in very good 
agreement with the well documented experimental value. 

Then, we examine the evolution of the exciton fine 
structure as a function of the exciton wavevector Q. Fig- 
ure[3]shows the calculated dispersion curves. For large Q, 
when the heavy hole light hole splitting becomes larger 
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than the exciton binding energy, the exciton splits into a 
"heavy" exciton formed of two twofold degenerate states 
and a "light" exciton formed of one threefold degener- 
ate plus one singlet states. Our calculation shows how 
energy levels interpolate between the small and large Q 
regimes. Finally, we note that when Q is along the [110] 
direction, our results show the full details of exciton state 
spin splittings, including both contributions of electron 
and hole spin splittings. 

In conclusion, we have devised a method that al- 
lows self-consistent definition of local wavefunctions 
within the extended basis empirical tight-binding model, 
and successfully used bulk exciton fine-structure as a 
parameter-free testbed. Extension to nanostructures is 
straightforward as long as bulk screening parameters 
are, like other tight binding parameters, transferable to 
nanostructures. This approach opens a route towards 
reconciling tight-binding and predictive evaluation of in- 
teractions between quasi-particles. 
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